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It is shown that single-peaked collisionless current sheets in a Harris-type equilibrium sponta- 
neously bifurcate as a result of chaotic scattering of electrons at fluctuating magnetic fields near the 
center of the sheet, as demonstrated by a 2D kinetic particle- in-cell simulation. For this effect to be 
simulated explicit particle advancing is necessary, since the details of the electron motion have to be 
resolved. Unlike previous investigations of triggering bifurcated current sheet (BCS) where initial 
perturbations or external pressure was applied the bifurcation is spontaneous if thermal noise is 
taken into account. A spontaneous current sheet bifurcation develops quicker than a tearing mode 
or other plasma instabilities. It is shown that in the course of the current sheet bifurcation the 
Helmholtz free energy decreases while the entropy increases, i.e. the new, bifurcated current sheet 
is in a more propable state than the single-peaked one. 

PACS numbers: 05.70.-a, 05.70.Ce, 52.35.Ra: 



The stability and possible unstable decay of current 
sheets play a central role in astrophyics as well as in 
the laboratory, e.g. for magnetic energy dissipation and 
reconnection For the investigations of current sheet 
stability often single-peaked current sheets are used, as 
the one derived by Harris [2j. 

The free energy of equilibria can cause a number 
plasma instabilities which spontaneously arise from ther- 
mal noise like the dissipative tearing mode instability in 
resistive [3j and collisionless plasmas Q. In fusion plas- 
mas spontaneous magnetic reconnection is found recently 
in the poloidal current sheet, which is an initial equilibria 
in the Reverse Field Pinch (RFP)[5[- Since the single- 
peaked current sheet equilibria is unstable and magnetic 
reconnection frequently takes place, the equilibria stabil- 
ity is an important issue for current sheet evolution. 

Magnetosphere is a natural plasma laboratory for the 
evolution of current sheet equilibrium. Indeed, more re- 
cent detailed investigations of current sheets, has shown 
that current sheets frequently are bifurcated (BCS) in- 
stead of sine; le-peaked [6] [7] [8]. It was also found that 
these BCS are electron dominated, e.g. by statistical 
analyses of measurements onboard the CLUSTER spac- 
craft mission [9]. 

BCS were also found in numerical simulations. In mag- 
netic reconnection plan BCS was interpreted, e.g., as a 
pair of slow mode shocks which develop in the recon- 
nection outflow region [To|[II|. However, BCS were also 
observed in minimum plasma inflow conditions, for which 
magnetic reconnection is not expected [12]. In the cur- 
rent direction, BCS formation without plasma inflow is 
discovered as a result of anomalous momentum transport 
due to pressure-gradient driven lower hybrid drift insta- 
bility [ll. 

In the direction perpendicular to current drift, current 
sheet splitting similar to BCS is observed after the sat- 
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uration of tearing mode instability [14]. Note that these 
authors used an implicit numerical scheme for advanc- 
ing the particles for their 2D particle-in-cell (PIC) code 
simulations. In order to initialize the tearing mode insta- 
bility quickly they also imposed perturbations to trigger 
instability growth. 

Schindler and Hesse [HI performed a one-dimensional 
particle-in-cell (ID PIC) simulation with an initial 
boundary pressing. They concluded that a quasisteady 
boundary compression forces a single-peaked current 
sheet evolves toward BCS, as an equilibrium relaxation 
process but not due to plasma instability 

We now have found that single-peaked collisionless 
current sheets might spontaneously bifurcate, without 
boundary compression or imposing tearing mode pertur- 
bation initially, as long as the electron thermal fluctua- 
tions are considered properly. We found that magnetic 
field fluctuations can initially start from thermal noise 
at the center of current. The magnetic field fluctuations 
may lead to chaotic scattering of the electrons out of the 
central (current peak) region of the sheet [16] . reducing 
the electron current flow there and adding current flows 
away from the center of the sheet as already discussed 
for laminar current sheets fl7j . 

To prove this hypothesis quantitatively we carried out 
2D electromagnetic particle-in-cell code simulations to 
investigate the evolution of Harris current sheet equi- 
librium. An explicit numerical scheme (XOOPIC) was 
implemented since the details of the electron cyclotron 
trajectories have to be calculated properly. The fastest 
processes up to electron plasma time scales uj~} = 

\J € m e /ne 2 is resolved. The chosen current sheet has 
a width 2 A = l.lbdi which is slightly larger than the 
ion inertial length di = cj uj p i to cover the full width 
of the ion dissipation region of the current sheet. The 
equilibrium magnetic field of a Harris sheet equlibrium is 
B y (x) = Botanh(x/ A) = ^/Afiks T Notanh(x / A) , where 
hs, T = Ti = T e , £?o and No are the Boltzmann con- 
stant, ion/electron temperatures, asymptotic magnetic 
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FIG. 1. The 2/- integrated current profiles (left panel) and current density at four different simulation times between t = and 
t = 8371, 2 cj^e 1 . The upper right panels show the total current and lower panels electron current densities. 



field and number density at the center of the sheet, re- 
spectively. The ratio of electron plasma frequency to 
electron cyclotron frequency is uj pe / uj ce = 2.87. The ion 
to electron mass ratio used is rrii/m e = 180, which is suf- 
ficient to separate the ion/electron motions. The grid size 
is of the Debye length dx = dy = X De = (eoT/Ne 2 ) 1 ^ 2 , 
which has sufficient spatial resolution of the detailed elec- 
tron motion. The simulation domain in the direction per- 
pendicular to the current sheet was L x — 13.3g^ and in 
the current sheet direction L y = 26.6^. This choice of L y 
allows a free development of the fastest growing tearing- 
type and other eigenmodes of the sheet (see, e.g., Eq.(23) 
in [HI). The boundary conditions for the particles and 
fields were chosen to be periodic in the y direction and 
conducting walls with particles reflection in the x direc- 
tion. 

The time step used was dt = 0.0872 uj~^ = 0.03 uo~^ to 
fulfill Courant condition and to resolve the detailed elec- 
tron oscillation and cyclotron motions. Note that this 
time step is much smaller than the one used in implicit 
numerical schemes, where, e.g., in [lij dt = 0.1 uj~1 ~ 
1.34 ^0.5 a;" 1 which does not track down the elec- 
tron oscillation and cyclotron motion. The later is signif- 
icant for magnetic scattering, which we will discuss soon. 
No initial perturbation was imposed as in [14]. Also, no 
boundary compression was imposed as in the ID PIC 
simulation of BCS formation (see Fig.l in jl5|). 

The left panels of FigJT]show the y integrated total cur- 
rent density I z (x) = I z , e {x) + I z ,i(x) and the individual 
contributions of electron I z e and ion I z i : for the ini- 
tial, single-peaked Harris-equilibrium currents by a thin 
solid line, for t = 4883.2 uj~} by a dashed line and for 
t = 8371.2 uj~c by a thick solid line. At the late stage 
t = 8371.2 uj~e ~ 16.212 cj" 1 two current peaks are 
formed which are well separated from each other. As one 
can see already in the middle panel of the left FigJT]the 
reduction of I z (x) near the center of the current sheet 
is mainly due the spatial redistribution of the electron 
current. Since the ion current is not changing much 
before t — 8371.2 uj~^ the right panels of Fig JT] depict 
the evolution of the total current density distribution in 
the x,y plane with time in the upper row and the re- 
sponsible for the current redistribution electron part in 




FIG. 2. Absolute values of the fluctuating \B X \ magnetic fields 
at the center of the current sheet, normalized to the asymp- 
totic magnetic field of the Harris equilibrium (left panel) and 
the thermal electron n{y) values at center of the current sheet 
(right panel) for the same four simulation times as in Figure 
1, right panels. 



the lower panels. While initially both the total and the 
electron currents are concentrated near the center of the 
sheet (first column in the right part of FigH]), already at 
t = 4883.2 uo~^ = 9.5 uj~^ = 1.5r C i a clear dip of the 
current density has developed around the center of the 
sheet. Here r C i = 2tt uj~^ denotes a cyclotron period. 
Hence, the bifurcation takes place at the time scale of 
an ion cyclotron period. This result indicates that the 
ion-to-electron mass ratio rrii/m e = 180 used in the sim- 
ulation the electron and ion dynamics are sufficiently well 
separated. 

Note further that the bifurcation due to boundary 
pressing takes place only after t ~ 400 uj~^ [l5[. Their 
mass ratio is rrii/m e = 25 and the bifurcation time in 
electron plasma period is 400 cj" 1 w 4.4 x 10 4 uj~J-. 

In order to understand this phenomenon one has to 
realize that plasmas are coarse-grained by their particles 
whose thermal motion lets currents and magnetic fields 
fluctuate at small scales. Hence, while the self-consistent 
magnetic field B y of a Harris-equilibrium vanishes at the 
center of the sheet, fluctuating magnetic fields remain. 

The left panel of Fig |2] shows \B x (x = 0)|, the absolute 
value of the normal field B x component of the magnetic 
field at the center of the sheet, divided by the asymptotic 
Harris sheet field B yo = Bq for the same simulation times 
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as chosen for the right panels in Fig. 1. As one can 
see, \B x (x = 0)| is finite most of the time unless the 
particles' thermal velocity is not taken into account as 
in the the uppermost column (t = 0). It is well known 
that particles transiting the center of current sheets can 
be chaotically scattered of the curvature of the magnetic 
field is comparable to the Larmor radii. According to 
Buchner and Zelenyi [16| the scattering is strongest, when 
the parameter k = \J RminPmax is unity, where R m in and 
rho max are the minimum curvature of the magnetic field 
and and the maximum particle gyroradius at the current 
sheet center. The right panel of Fig|2]depicts the n values 
for thermal electrons along the y direction for the same 
moments of time for which the B x (y) fields are shown 
in the left panel of this Figure and in the right panel of 
Fig.l. As one can see in the right panel of Fig|2]before t — 
4883 k reaches unity only at a few positions, i.e. only 
at a few places the electrons are strongly scattered. After 
t = 4883 electrons are strongly scattered chaotically 
at many positions. 

The consequences of this strong scattering can be de- 
scribed best by means of the action integral of the fast 
motion motion / = §v z dz [19|. For n < 1 this ac- 
tion integral can be expressed in its normalized form 

as V = ( 2 fc 2 -i ) /a/bWj where k(y') is a function 
of y' , the appropriately normalized y coordinate of the 
slowly moving guiding center and are two dif- 

ferent functions of complete elliptic integrals for k < 1 
and k > 1, respectively [16|. For n < 1 the action in- 
tegrals V are adiabatically conserved, i.e. integrals of 
motion, as long as they stay away from the separatrix 
which is reached when k —> 1. For k < 1 such quasi- 
adiabatic orbits cross the center of the current sheets 
while for k > 1 they will gyrate at some distance from 
the sheet center. Particles with 1 < V < I' max = 1.16 
are trapped on crossing the current sheet orbits k < 1. 
But in thin current sheets most particles will have I' < 1, 
which enables them to eventually reach a separatrix in 
the velocity space, i.e. change from the region k < 1 
to k > 1, at k = 1 changing from meandering across 
the sheet, drifting in the direction of the original sheet 
current to a gyration away from the current sheet center 
where they drift in the opposite direction, causing dia- 
magnetic currents away from the sheet center [16| . While 
for n <C 1 particles the value of the quasi-integral of mo- 
tion V stays practically unchanged during the sepratrix 
encounter, for n — > / it can be essentially changed. If the 
obtained value of V is very small only a small amount of 
kinetic energy is left in the perpendicular to the magnetic 
field velocity direction while most energy is in the slow 
drift motion. These particles contribute significantly to 
the built up of diamagnetic currents away from the sheet 
center. For small V the asymptotic expressions for the 
elliptic integrals reveal V « 3/1 671-/^ ft 6 , where k tp cor- 
respond to the turning point of the drift motion in y' 
(v y t = 0). From the condition for turning the drift in 
the y direction one obtains the position of the turning 
point by solving the equation n 2 y[ p = 2k 2 p — 1. Solving 
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FIG. 3. Profiles of electron drift velocity V e , z (upper panel) 
and number density N e (lower panel), both integrated along 
the y axis. 



this equation for k tp s and calculating the mean z-position 
for the turning point location one obtains a the distance 
Az\J\p t he{^§K /3) 2 I' the from the current sheet at which 
most of the particles drift in the dia-magnetic current 
direction. Here I' the means the quasi-adiabatic invariant 
for typical thermal electrons, the bulk drift velocity of 
the electrons is much smaller and can be neglected. With 
the time being Az will become the position of the maxi- 
mum of the dia-magnetic flow of the electrons where the 
initial current profile is modified most by electron flows 
while the electron flow at the center of the current sheet 
is reduced. For the parameters used in the simulation 
one obtains Az = 27\De- This theoretically predicted 
distance corresponds to the one found in simulation. 

The Y integrated electron drift V e , z and density N e 
profiles are show in Fig J3] for three moments of time. Note 
that such current profile was obtained by Hoshino et al. 
[7j from a statistical analysis of many current sheet en- 
counters in the Earth's magnetotail. The electron num- 
ber density does not change much in the course of the 
bifurcation, which happens, instead, in the electron drift 
velocity space. This is consistent with the finding that 
the current depletion at the center is caused mainly by 
the redistribution of the electrons by pitch-angle scat- 
tering in the velocity space when crossing the sepratrix 
between meandering (k < 1) and gyrating away from the 
sheet center (k > 1). Due to the mass ratio only the 
electrons undergo strong chaotic scattering thin Harris- 
type current sheets. The new, bifurcated current sheet is 
again in equilibrium. In fact there is an infinite number of 
possible current sheet equilibria. A number of analytical 
[20[ and non-analytical equilibrium solutions have been 
found which are more realistic for space current sheets 
than Harris sheets [2l|. For a survey see, e.g., [22|. A 
current sheet bifurcated out of a single peaked Harris cur- 
rent sheet by the chaotic electron scattering is close to 
the equilibrium found by Camporeale and Lapenta [lij , 
their case (b). Note that the BCS, naturally obtained 
via chaotic electron scattering, is more probable than the 
single peaked Harris sheet. This can be demonstrated by 
calculating the entropy of the system. For a plasma with 
a continuous particle distribution, it is appropriate to 
consider the relative Kullback-Leibler entropy which is 
always positive [23|. With respect to a reference distri- 
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FIG. 4. Evolution of the total relative entropy (left panel) and 
of the electron (upper right panel) and ion entropies (lower 
right panel) separately. 

bution q(y) at t = it can be written as 

S KL (t)= dvf(v,t)ln(^^) 

The relative entropy as a sum of the electron and ion 
contributions grows increases in the course of current 
sheet bifurcation. Looking at the electrons and ions sepa- 
rately one can see that the entropy of electrons is, indeed, 
steadily increasing after about t = 3800 uSl} while the ion 
entropy oscillates (right columns of Fig|4j). The stability 
of a steady state equilibrium can be analyzed by calcu- 
lating the Helmholtz free energy F = U — T S [24[ . Here 
the U is the internal energy, a sum of the particles kinetic 
energy and the field energy, T is the temperature energy 
and S is the entropy. Less free energy F corresponds to 
a more stable equilibrium. In the course of the bifurca- 
tion, the internal energy U is conserved while T and S are 
increasing. Hence the bifurcated current sheet contains 



less free energy. This is why it is more favored in nature 
than single-peaked equilibria like the Harris current sheet 
and it is more stable. 

2D PIC simulations confirmed that an initially sin- 
gle peaked Harris equilibrium current sheet does spon- 
taneously and quickly bifurcates a natural consequence 
of chaotic electron scattering due to their fluctuations at 
the center of the sheet. The bifurcation is faster than 
spontaneous plasma instabilities as of the tearing mode. 
The single peaked sheet bifurcates without initial per- 
turbations or imposed external pressure. The bifurcated 
state is more favorable because electrons gaining after 
chaotic scattering at the sheet center toward small val- 
ues of the quasi-adiabatic invariant of motion V spends 
longer time in the cucumber phase of their motion away 
from the sheet center rather than meandering at the sheet 
center. Therefore it can be simulated only by numerical 
schemes that explicitly resolve the details of the electron 
motion. The spontaneous bifurcation of single peaked 
current sheets explains the frequent observation of bifur- 
cated current sheets in quite situations without plasma 
inflows and reconnection as they are more stable than 
single peaked sheets. Their influence on plasma insta- 
bilities is stabilizing, delaying, lowering the growth rate, 
e.g., of the tearing mode instability. 
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